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Abstract 

In  this  paper,  we  propose  a  novel  framework,  called  Dinkelbach  NCUT  (DNCUT),  which 
efficiently  solves  the  normalized  graph  cut  (NCUT)  problem  under  general,  convex  constraints, 
as  well  as,  under  given  priors  on  the  nodes  of  the  graph.  Current  NCUT  methods  use  generalized 
eigen-decomposition,  which  poses  computational  issues  especially  for  large  graphs,  and  can  only 
handle  linear  equality  constraints.  By  using  an  augmented  graph  and  the  iterative  Dinkelbach 
method  for  fractional  programming  (FP),  we  formulate  the  DNCUT  framework  to  efficiently 
solve  the  NCUT  problem  under  general  convex  constraints  and  given  data  priors.  In  this  frame¬ 
work,  the  initial  problem  is  converted  into  a  sequence  of  simpler  sub-problems  (i.e.  convex, 
quadratic  programs  (QP’s)  subject  to  convex  constraints).  The  complexity  of  finding  a  global 
solution  for  each  sub-problem  depends  on  the  complexity  of  the  constraints,  the  convexity  of  the 
cost  function,  and  the  chosen  initialization.  However,  we  derive  an  initialization,  which  guaran¬ 
tees  that  each  sub-problem  is  a  convex  QP  that  can  be  solved  by  available  convex  programming 
techniques.  We  apply  this  framework  to  the  special  case  of  linear  constraints,  where  the  solution 
is  obtained  by  solving  a  sequence  of  sparse  linear  systems  using  the  conjugate  gradient  method. 
We  validate  DNCUT  by  performing  binary  segmentation  on  real  images  both  with  and  without 
linear/nonlinear  constraints,  as  well  as,  multi-class  segmentation.  When  possible,  we  compare 
DNCUT  to  other  NCUT  methods,  in  terms  of  segmentation  performance  and  computational  ef¬ 
ficiency.  Even  though  the  new  formulation  is  applied  to  the  problem  of  spectral  graph-based, 
low-level  image  segmentation,  it  can  be  directly  applied  to  other  applications  (e.g.  clustering). 


1  Introduction 

This  paper  is  about  extending  normalized  graph  cuts  so  the  result  satisfies  general,  convex  constraints 
under  given  priors  on  the  graph.  The  new  formulation  is  applied  to  and  presented  in  terms  of  low 
level  image  segmentation  using  spectral  graph  theory,  a  problem  that  has  received  extensive  attention 
recently.  However,  this  framework  can  be  applied  to  general  data  clustering  problems,  where  the 
aforementioned  segmentation  problem  is  a  special  case.  In  this  problem,  an  image  is  represented  by 
an  undirected  graph,  wherein  nodes  correspond  to  image  pixels  or  regions,  and  edges  connect  pairs 
of  nodes.  An  edge  has  a  weight  proportional  to  the  similarity  of  the  properties  of  the  connected  nodes 
(e.g.  pixel  intensities).  Given  such  a  graph  representation,  image  segmentation  becomes  equivalent 
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to  partitioning  the  nodes  of  the  graph  into  disjoint  sets,  or  segments,  which  optimize  a  given  cost 
function.  Such  a  partition  is  denoted  as  a  graph  cut.  Traditionally,  the  cost  of  a  cut  between  two 
segments  is  the  sum  of  weights  corresponding  to  the  graph  edges  that  need  to  be  severed  to  produce 
this  segmentation.  However,  other  graph  cut  formulations  exist,  so  graph  cut  methods  are  categorized 
by  the  global  cost  function  they  optimize.  If  the  objective  is  to  minimize  the  cost  of  a  cut  between 
two  or  more  segments,  then  the  problem  can  be  formulated  as  a  min-cut  and  efficiently  solved  as 
a  max-flow  problem  [3, 5, 11, 18].  This  cost  can  include  node  priors  by  introducing  artificial  nodes 
corresponding  to  the  number  of  desired  segments.  The  edge  weights  from  these  nodes  to  the  rest  of 
the  graph  embed  the  node  priors.  In  what  follows,  we  will  use  the  label  unnormalized  graph  cuts  to 
denote  this  graph  cut  problem. 

Despite  the  merits  of  the  unnormalized  graph  cut  formulation  and  its  proposed  algorithms,  it  is 
biased  towards  producing  cuts  that  contain  a  small  number  of  nodes.  Consequently,  a  normalized 
version  of  this  cost  function  has  been  proposed,  in  which  the  cost  of  the  normalized  cut  (NCUT)  in¬ 
corporates  the  association  cost  of  each  segment  (i.e.  the  sum  of  edge  weights  connecting  nodes  of  this 
segment  to  all  other  nodes  in  the  graph),  thus,  suppressing  the  appearance  of  small  cuts  [23].  For  a 
detailed  comparison  of  these  two  formulations,  refer  to  [25].  Minimizing  the  resulting  NCUT  objec¬ 
tive  function  is  an  NP  hard  discrete  optimization  problem,  so  it  is  relaxed  to  take  on  real  values.  The 
popularity  of  this  method  can  be  partially  attributed  to  its  closed  form  approximation  of  the  optimal 
solution,  using  generalized  eigen-decomposition.  Unconstrained  NCUT  has  been  used  extensively 
for  image  segmentation  and  data  clustering.  In  [26],  the  authors  constrain  the  NCUT  problem  by 
incorporating  additional  grouping  constraints,  in  the  form  of  linear  homogeneous  equalities.  This  is 
extended  to  the  case  of  non-homogeneous  equalities  in  [9]. 

In  all  the  above  NCUT  methods  [9, 23, 26],  the  NCUT  solution  is  obtained  from  either  a  single 
generalized  eigen-decomposition  or  a  sequence  of  such  decompositions.  Even  though  these  methods 
offer  certain  advantages,  their  following  shortcomings  need  to  be  addressed.  (1)  These  methods  only 
allow  for  linear  equality  constraints.  This  is  due  to  their  underlying  use  of  eigen  decomposition  to 
minimize  the  Rayleigh  quotient.  For  example,  [6]  showed  that  it  is  NP-hard,  in  general,  to  solve  for 
eigenvectors  under  linear  inequalities.  This  suggests  that  if  a  unifying  framework  for  such  constrained 
NCUT  problems  is  desired,  it  should  avoid  using  this  quotient  formulation.  (2)  Unlike  unnormalized 
graph  cut  techniques,  these  methods  do  not  embed  any  unary  terms  (i.e.  priors  on  individual  graph 
nodes)  in  the  NCUT  cost  function,  which  is  equivalent  to  assuming  a  zero  prior  for  all  nodes.  (3) 
Since  these  methods  compute  generalized  eigenvectors  of  large  matrices  and  form  null  spaces  of 
highly  rank  deficient  matrices,  their  computational  complexity  remains  a  practical  issue,  despite  the 
special  measures  that  were  considered  for  complexity  reduction. 

Contributions 

In  this  paper,  we  present  a  first  step  in  the  direction  of  formulating  a  unifying  framework  for  convexly 
constrained  NCUT  problems,  which  addresses  the  aforementioned  limitations.  The  contributions  of 
this  framework  are  threefold,  (i)  It  allows  the  efficient  solution  of  the  NCUT  problem  under  general 
convex  constraints.  It  uses  the  Dinkelbach  method  to  transform  the  initial  fractional  problem  into  a 
sequence  of  convexly  constrained,  quadratic  programming  (QP)  problems,  whose  convexity  is  en¬ 
sured  by  a  suitable  initialization  that  we  construct.  More  importantly,  the  global  solutions  of  these 
problems  converge  superlinearly  to  the  required  solution  of  the  fractional  problem.  Since  the  Dinkel¬ 
bach  method  is  central  to  our  framework,  we  denote  our  proposed  method  as  Dinkelbach  NCUT 
(DNCUT).  In  fact,  the  Dinkelbach  method  was  recently  used  in  the  context  of  parametric  max-flow 
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algorithms  in  [12].  However,  the  energy  function  to  be  minimized  was  unconstrained  and  hence  no 
effort  was  made  to  construct  a  valid  initialization  that  allows  efficient  estimation  of  the  global  so¬ 
lution  for  each  Dinkelbach  iteration,  (ii)  Our  framework  can  incorporate  prior  knowledge  of  nodes 
belonging  to  different  segments,  which  is  equivalent  to  the  unary  term  present  in  the  energy  function 
minimized  by  unnormalized  graph  cut  algorithms.  This  case  arises  in  various  supervised  problems 
including  interactive  segmentation  and  supervised  clustering.  These  problems  can  not  be  solved  us¬ 
ing  the  current  NCUT  algorithms,  (iii)  Guaranteeing  the  convexity  of  each  Dinkelbach  sub-problem 
allows  this  framework  to  handle  general  convex  constraints,  due  to  the  availability  of  efficient  con¬ 
vex  optimization  techniques.  Therefore,  our  proposed  algorithm  refrains  from  eigen-decomposition 
and,  in  the  case  of  linear  constraints,  it  degenerates  into  solving  a  set  of  sparse  linear  systems  using 
conjugate  gradients. 

Mathematical  Notation 

Before  proceeding  with  the  details  of  the  paper,  we  summarize  the  mathematical  notations  used  in 
the  remainder  of  the  paper.  In  what  follows,  a  matrix  is  denoted  by  a  bold,  uppercase  letter  (e.g.  W). 
A  vector  is  denoted  by  a  bold,  lowercase  letter  with  an  arrow  above  it  (e.g.  x).  A  scalar  is  denoted 
by  a  normal  lowercase  letter.  Elements  of  vectors  and  matrices  are  indexed  using  parentheses.  For 
example,  x(i)  represents  the  ith  element  of  x,  while  W (i,j)  represents  the  element  of  W  at  the  ith 
row  and  the  jth  column.  Furthermore,  a  parenthesized  superscript  is  added  to  a  variable  to  denote  its 
value  at  a  given  iteration.  For  example,  represents  the  value  of  x  at  the  kth  iteration.  Also,  every 
optimization  problem  has  a  solution  appended  with  a  superscript  *  (e.g.  x*)  and  constraints  (if  any) 
listed  underneath  the  cost  function  prefixed  by  s.t.,  the  abbreviation  of  such  that. 

We  separate  the  discussion  of  applying  hard  and  other  general  convex  constraints  to  the  NCUT 
problem.  This  is  done  because  hard  constraints  have  direct  impact  on  the  cost  function  itself,  while 
the  other  constraints  do  not.  In  Section  2,  we  consider  the  problem  of  applying  hard  constraints  (e.g. 
placed  by  the  user  via  interactive  placement  of  constraints,  which  we  call  here  interactive  NCUT). 
In  Section  3,  we  describe  the  augmented  graph  structure  used  in  formulating  and  solving  the  NCUT 
problem.  In  Section  4,  we  use  the  Dinkelbach  approach  [8]  to  formulate  a  unifying  framework  for 
solving  the  convexly  constrained  NCUT  problem.  We  also  describe  an  algorithm  to  solve  this  problem 
with  linear  constraints.  Finally,  we  test  the  algorithm  by  applying  it  to  low-level  segmentation  of  real 
images  in  Section  5. 

2  Hard  Constraints 

In  this  section,  we  consider  the  NCUT  problem  under  the  first  type  of  convex  constraints,  namely 
hard  constraints.  To  put  this  problem  into  context,  we  visit  the  unnormalized  graph  cut  problem, 
which  will  also  help  us  introduce  the  interactive  version  of  the  NCUT  problem.  Due  to  the  nature 
of  the  two  different  cost  functions  optimized  in  the  unnormalized  graph  cut  problem  and  its  corre¬ 
sponding  NCUT  problem,  the  analysis/solution  of  these  two  problems  is  quite  different.  However, 
we  still  make  use  of  two  general  concepts  used  in  the  unnormalized  graph  cut  context:  the  augmented 
graph  structure  and  the  notion  that  hard  constraints  directly  affect  the  graph’s  weight  matrix.  We  will 
elaborate  on  these  two  points  in  what  follows.  Next,  we  introduce  the  terminology  to  be  used  and  give 
a  brief  description  of  the  unnormalized  graph  cut  and  interactive  unnormalized  graph  cut  problems. 
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Figure  1 :  Example  of  a  graph  cut  (x),  which  separates  nodes  1,2,3  (labeled  as  A)  from  4,  5,  6  (labeled 
as  B).  The  cost  of  x  is  the  sum  of  the  edge  weights  that  have  to  be  cut  to  produce  the  final  labeling. 
Note  that  the  regional  costs  are  represented  by  the  perforated  blue  lines,  while  the  solid  green  lines 
represent  the  boundary  costs. 

2.1  Unnormalized  Graph  Cuts  and  Hard  Constraints 

Let  Q  =  {V,  £)  be  an  undirected  graph,  where  V  denotes  the  set  of  nodes  contained  in  Q,  and  £ 
the  set  of  corresponding  edges.  Let  W  denote  the  edge  weight  matrix  of  Q,  where  W  (z,  j)  is  a 
non-negative  weight  associated  with  the  edge  connecting  nodes  i  and  j  in  V.  This  non-negativity  is 
the  only  restriction  on  W,  whereby  the  computation  of  pairwise  edge  weight  values  is  application- 
dependent.  When  applied  to  low-level  spectral  graph  image  segmentation,  an  edge  weight  between 
two  pixel  nodes  is  commonly  based  on  the  similarity  in  color  and/or  image  gradient  between  these 
two  pixels.  A  graph  cut  x  is  a  binary  segmentation  of  the  nodes  in  V ,  where  Vfc  =  1,  •  •  •  ,  \V\, 
x(fc)  =  1  if  the  kth  node  belongs  to  segment  A  and  x(fc)  =  0  if  it  belongs  to  segment  B  (refer 
to  Ligure  1).  The  aim  of  a  graph  cuts  algorithm  is  to  find  x  that  minimizes  a  given  cost  function. 
In  [11],  the  cost  function  was  defined  as  E^(A^B)  =  \R(x)  +  S'(x).  The  regional  cost  term, 
#(*)  =  Eie-P  Ri(x(i)),  is  the  sum  of  costs  incurred  by  assigning  each  node  i  to  its  label  x(z).  The 
boundary  (cut)  term,  S'(x)  =  cut(A,B)  =  J2ieA  jeB^  (hj)>  the  cost  of  the  cut  resulting  in 
segments  A  and  B.  Numerous  efficient  algorithms  have  been  developed  to  find  the  global  minimizer 
of  Ez(A,B)  [3]. 

The  interactive  unnormalized  graph  cuts  problem  is  an  extension  of  the  aforementioned  prob¬ 
lem  with  the  addition  that  some  hard  constraints  are  applied  to  the  nodes  of  V  (i.e.  labels  of  some 
nodes  are  known  beforehand).  These  constraints  can  originate  from  direct  user  interaction  or  domain 
specific  knowledge.  Many  recent  works  have  addressed  this  problem  under  different  forms  of  user 
interaction  [4, 14,22].  We  denote  Sa  and  Sb  to  be  the  sets  of  nodes  satisfying  the  hard  constraints 
(i.e.  whose  labels  are  known  beforehand  to  be  1  or  0,  respectively).  Under  these  constraints,  the 
minimization  problem  becomes  more  difficult  to  solve.  However,  Boykov  et.  al.  [4]  proved  that  there 
is  no  need  to  explicitly  solve  this  new  problem.  They  show  that  it  is  equivalent  to  an  unconstrained 
unnormalized  graph  cut  problem  on  Q  with  the  edge  weights  appropriately  modified  to  implicitly 
reflect  the  hard  constraints.  While  such  an  equivalence  exists  for  the  interactive  unnormalized  graph 
cut  problem,  it  does  not  transfer  to  the  corresponding  NCUT  problem,  which  we  will  appropriately 
call  the  interactive  NCUT  problem.  Next,  we  show  how  these  hard  constraints  are  explicitly  applied 
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to  this  NCUT  problem.  We  study  this  case  independently,  since  it  will  have  direct  influence  on  the 
cost  function  to  be  minimized. 


2.2  NCUT  and  Hard  Constraints 


A  normalized  cut  of  Q  into  segments  A  and  B  has  the  following  cost:  NCUT^(A,  B) 


cut(A,B) 


 cut(A,B) 


jieAjev 


ZieBjevW^j) 

formulation  defined  in  [23]  is: 


W  (i,j) 

,  where  the  cut  cost  is  normalized  by  the  association  cost  of  each  segment.  The  NCUT 


y NCUT  =  arg  min 


yT(D  -  W)y 

yTDy 


s.t. 


y (i)  e  {1  ,-bj  V*  =  1,  •  • 
yTDl  =  0 


\v\ 


(1) 


where  D  is  the  degree  matrix  of  Q  (i.e.  D  =  diag  (\'Vi'j )  and  b  =  ^x(,)>0p|]’^.  Here,  we  note 

that  b  quantifies  how  connected  (similar)  each  segment  is  to  the  rest  of  the  graph.  A  large  value  of  b 
(>  1)  means  that  the  nodes  of  segment  A  have  stronger  connections  to  the  rest  of  the  graph  than  those 
of  segment  B.  Since  the  value  of  b  is  dependent  on  the  final  labeling,  solving  this  problem  exactly 
becomes  infeasible.  In  general,  the  optimization  problem  in  Eq  (1)  is  NP  hard,  so  it  is  relaxed  to 
render  a  real  solution.  This  vector  is  discretised  as  a  post  processing  step,  which  does  not  incorporate 
the  value  of  b.  In  [23],  the  authors  show  that  the  solution  to  the  relaxed  problem  can  be  obtained 
in  closed  form  by  solving  a  generalized  eigenvalue  problem  (i.e.  eig(D  —  W,  D)).  This  is  a  direct 
conclusion  from  the  fact  that  this  relaxed  problem  is  in  the  form  of  a  Rayleigh  quotient.  In  fact, 
Yncut  is  computed  as  the  generalized  eigenvector  corresponding  to  the  second  smallest  eigenvalue. 
In  [23],  the  NCUT  framework  was  applied  to  low-level  image  segmentation  and  the  edge  weight 
between  each  pair  of  nodes/pixels  was  computed  using  intervening  contours  [16]. 

Now,  let  us  extend  Eq  (1)  to  include  some  hard  constraints,  rendering  the  interactive  NCUT  prob¬ 
lem.  To  the  best  of  our  knowledge,  the  interactive  NCUT  problem  has  not  been  addressed  previously 
in  the  literature.  To  solve  this  problem,  we  cannot  use  the  same  graph  used  in  the  traditional  NCUT 
formulation.  This  graph  is  composed  solely  of  pixel  nodes  originating  from  the  image  itself.  To  this 
graph,  we  add  two  artificial  nodes  A  and  B  to  represent  the  two  segments,  thus,  differentiating  them 
from  the  pixel  nodes.  This  augmented  graph  will  allow  us  to  incorporate  hard  constraints  and  in¬ 
clude  prior  knowledge.  We  decompose  the  labeling  vector  y  into  three  disjoint  parts:  (1)  y+  (of  size 
N+  =  \Sa\),  which  corresponds  to  the  pre-labeled  nodes  of  Sa,  (2)  y_  (of  size  =  \Sb |),  which 
corresponds  to  the  pre-labeled  nodes  of  Sb,  and  (3)  yu  (of  size  Nu),  which  designates  the  unknown 
labels  of  the  rest  of  the  nodes  in  the  graph.  We  use  the  augmented  graph  in  Figure  2  to  illustrate  an 
example.  We  update  W  and  D  to  include  the  hard  labels,  as  shown  below.  The  nodes  corresponding 
to  y+  are  listed  first,  followed  by  those  of  y_,  and  then  yu.  Similar  to  the  interactive  unnormal¬ 
ized  graph  cut  problem  seen  before,  we  embed  the  hard  constraints  into  the  edge  weights  by  setting 
W(i,  j)  =  W  (V,/)  =  K,  W  (*,*')  =  W  (j,/)  =  0  for  every  node  pair  ( i,j )  g  Sa  x  Sa 

and  (/, /)  €  Sb  x  Sb-  This  means  that  nodes  belonging  to  the  same  segment  are  maximally  sim- 
ilar,  where  the  similarity  value  is  a  constant  denoted  by  K.  For  now,  K  and  b  are  left  unknown.  In 
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Figure  2:  Example  of  a  graph  with  hard  constraints  colored  in  red.  Here,  N+  =  N_  =  2,  where  two 
pixel  nodes  are  already  labeled  as  A  and  B.  The  perforated  blue  edges  are  designated  weights  that 
embed  prior  knowledge  of  the  unknown  nodes.  The  bold  red  edges  have  weights  of  K.  Bold  green 
edges  represent  the  similarity  between  the  unknown  nodes  they  connect. 


Section  4.2.2,  we  show  how  they  are  computed.  Similarly,  the  nodes  belonging  to  different  segments 
are  minimally  similar.  Furthermore,  we  set  W  (i,j)  =  p+  (i)  and  W  (hj'j  =  Pu  W  f°r  everY 

i  0  Sa  U  Sb,  j  E  Sa,  and  j  G  Sb-  Here,  prior  knowledge  can  be  incorporated.  The  edge  weight 
between  an  unknown  node  and  the  hard  constraints  (i.e.  the  pre-labeled  nodes)  can  be  viewed  as  the 
likelihood  of  that  unknown  node  belonging  to  either  segment.  So,  if  a  likelihood  node  model  (e.g. 
a  Gaussian  Mixture  Model  (GMM))  is  available,  it  is  used  to  evaluate  p+  and  p“.  In  our  experi¬ 
ments,  given  hard  constraints,  we  assume  the  Nearest  Neighbor  model,  so  we  set  p+  (i)  and  p~  ( i ) 
to  the  maximum  edge  weight  between  node  i  and  nodes  of  Sa  and  Sb  respectively.  In  the  absence 
of  prior  knowledge,  these  edge  weights  are  set  to  the  same  constant.  In  our  experiments  and  in  the 
absence  of  hard  constraints,  we  set  p+  (i)  =  p~  (i)  =  since  all  edge  weights  (based  on  intervening 
contours  [15])  take  values  in  [0, 1]  for  the  case  of  low-level  image  segmentation.  Consequently,  the 
updated  W  and  D  matrices  are: 


K 1 

o  p+  ' 

"  (KN+  +  1T^+)I  0 

0 

w  = 

0 

p+T 

L  u 

K 1  P^ 

p^t  wttJ 

;  d  = 

0  (KN.  +  1T^)1 

0  0 

0 

Du  +  Pd 

where  PD  =  diag  (iV+p+  +  iV_ pu  ) ,  P+  =  1  (p+ )T,  and  Pu  =  1  (pu  )T  (refer  to  Figure  2). 


3  DNCUT  Graph  Structure 

As  opposed  to  the  traditional  NCUT  formulation  that  only  uses  pixel  nodes  to  construct  the  graph, 
we  adopt  the  augmented  graph  structure  used  in  the  unnormalized  graph  cut  problem.  As  illustrated 
before  in  Figure  2,  the  augmented  graph  adds  artificial  nodes  to  the  original  graph.  These  added  nodes 
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correspond  to  the  “sink”  and  “source”  terminals  in  [18].  This  augmented  graph  structure  allows  for 
two  major  benefits. 

1.  It  inherently  incorporates  prior  knowledge.  We  designate  p+  and  p“  to  be  the  vectors  of  edge 
weights  connecting  the  unknown  nodes  to  the  nodes  in  segments  Sa  and  Sb  respectively.  In 
fact,  they  correspond  to  the  unary,  regional  cost  terms  included  in  the  unnormalized  graph  cut 
formulation.  They  can  also  be  viewed  as  the  likelihoods  that  the  unknown  nodes  belong  to 
each  of  the  two  segments.  As  such,  the  traditional  NCUT  problem  becomes  a  special  case 
of  the  DNCUT  framework,  when  p+  =  p“  —>  0  (i.e.  zero  priors)  and  K  — ►  0.  Figure  3 
shows  an  example  of  how  prior  knowledge  on  the  nodes  of  the  augmented  graph  can  improve 
segmentation/labelling  results.  Figure  3(b)  is  the  binary  DNCUT  segmentation,  if  uniform  prior 
knowledge  is  assumed  on  the  nodes  of  this  graph.  The  prior  knowledge  in  Figure  3(c)  depicts 
how  nodes  belonging  to  class  Sa  and  Sb  should  look  like  in  the  image.  Note  the  significant 
difference  in  segmentation  when  comparing  Figure  3(b)  to  Figure  3(d).  Traditional  NCUT 
methods  cannot  incorporate  prior  knowledge  on  the  nodes  of  the  graph. 


(a):  original  image  (b):  DNCUT  (c):  prior  (d):  DNCUT  with  prior 


Figure  3:  For  the  image  in  (a),  the  unconstrained  binary  DNCUT  segmentation  is  provided  in  (b). 
This  segmentation  is  obtained  when  all  unknown  nodes  have  equally  likely  prior  probabilities  of 
belonging  to  Sa  and  Sb  (i.e.  p+  (i)  =  p“  (i)  =  \).  However,  (d)  represents  the  DNCUT  solution 
when  nonuniform  prior  probabilities  are  used.  These  probabilities  are  based  on  how  similar  the 
unknown  nodes  are  to  the  nodes  delineated  by  the  green  (Sa)  and  red  (Sb)  strokes.  We  used  a  simple 
GMM  (two  Gaussians)  likelihood  model.  Note  how  the  incorporation  of  prior  knowledge  rendered 
the  binary  DNCUT  segmentation  more  meaningful. 

2.  It  supports  an  indirect  connection  between  every  pair  of  pixel  nodes  in  the  graph,  while  preserv¬ 
ing  the  sparsity  of  W.  In  the  case  of  image  segmentation  and  due  to  memory  restrictions,  W 
is  made  sparse  by  setting  the  edge  weights  between  far  pixels  to  0.  For  the  traditional  NCUT 
problem,  this  can  yield  pairs  of  nodes  that  have  high  similarity  but  are  neither  directly  nor  in¬ 
directly  connected  in  the  graph.  This  biases  the  segmentation  to  assign  such  nodes  to  different 
segments.  This  usually  occurs  due  to  occlusions.  On  the  other  hand,  the  DNCUT  framework 
ensures  an  indirect  connection  between  these  nodes,  via  the  nodes  in  Sa  and  Sb ,  thus,  alleviat¬ 
ing  the  previous  segmentation  bias.  In  Figure  4,  we  show  an  example  of  the  traditional  NCUT 
and  the  DNCUT  solutions  to  the  unconstrained  NCUT  problem  described  in  Section  2.2  [23]. 
Note  that  these  solutions  are  binary  solutions  obtained  by  discretising  the  real  solutions  to  the 
unconstrained  NCUT  problem. 
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(a):  original  image 


(b):  traditional  NCUT 


(c):  DNCUT 


Figure  4:  For  the  image  in  (a),  the  traditional  NCUT  solution  to  the  unconstrained  NCUT  problem  is 
illustrated  in  (b)  and  the  DNCUT  solution  in  (c).  Note  that  similar  pixel  nodes  (i.e.  sky)  are  assigned 
to  different  segments  in  (b)  and  to  the  same  segment  in  (c). 


4  DNCUT  Framework  under  Hard  &  Convex  Constraints 

In  this  section,  we  highlight  the  details  of  our  proposed  DNCUT  framework.  We  consider  the  con- 
vexly  constrained  NCUT  problem,  where  hard  and/or  other  convex  constraints  are  applied.  Given 
sets  Sa  and  Sb  and  the  updated  weight  matrix  W,  we  formally  define  the  problem  as  a  fractional 
quadratic  program  (FQP)  with  convex  constraints,  as  shown  in  Eq  (2).  Note  that  the  following  setup 
is  the  same  even  if  no  hard  constraints  exist.  For  this  special  case,  Sa  and  S b  only  contain  a  single 
node  each  (i.e.  N+  =N-= 1). 


where 


:  arg  mm 


fuQyu  +  mTyu  +  (a  -  c) 


y«  Ry  u  +  a 


y u(*)e{i,-&}  V«  =  l,  •  •  •  ,NU 
fu' R1  +  <7  =  0 
$i( Yu)  <0,  Vi  =  1,  —  ,  |I| 
Ai(fu)  =  0,  Vj  =  !,•■■  ,|E| 


(2) 


Q  —  (Du  —  Wu)  +  Pd  —  Lu  +  Pc,  Qe 
m  =  2  (bN-p~  -  iV+p+N 


a  =  N+  (. KN+  +  1T^+)  +  b2N-  (kN-  +  lTp" 


c  =  K  (N+)  +(bN-) 
R  =  +  P D-,  FI  G  E>tr 


' Nu 

q  =  N+  [^KN+  +  lTp+)  -  bN _  (j<N-  +  lrp 
k$i(x)  and  'Pj(x)  are  convex  Vi  =  1  •  •  •  |I|,  j  =  1  •  •  •  |E| 


?T„- 


(3) 


We  define  the  Laplacian  matrix  corresponding  to  the  unknown  nodes  as  Lu,  which  is  known 
to  belong  to  =  {Xe  rNuxNu  :  x  =  XT,  X  ^  0|  (i.e.  the  set  of  symmetric  positive  semi- 
definite  matrices  of  size  Nu  x  Nu)  [20].  For  image  segmentation,  Wu  and  Q  are  sparse,  in  general. 
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Q>i(yu)  and  ^j{ yu)  are  general  convex  constraints  that  can  be  linear  (e.g.  partial  groupings  [26])  or 
non-linear  (e.g.  upper  bounds  on  |  |yn|  |2)- 

We  keep  b  and  K  as  scalar  variables.  Consequently,  m,  a,  c,  and  q  are  variables  too.  Note  that  the 
value  of  b  depends  on  y  *  itself  and  that  the  problem  in  Eq  (2)  is  still  an  NP  hard  discrete  optimization 
problem,  so  we  propose  two  forms  of  relaxation:  (1)  we  relax  yu  to  take  on  real  values  and  (2)  we 
assume  b  takes  on  a  constant  value  bo.  In  Section  4.2.2,  we  show  how  K  and  bo  are  computed  in  the 
DNCUT  framework.  The  relaxed  problem  becomes  the  FQP  defined  in  Eq  (4). 


y  *  =  arg  min 


H(  y«) 


G(  yu)_ 


=  arg  min 


ylQyu  ±  mry u  +  {a -c) 

y«Ry«  +  a 


f  y^Ri  +  q  =  o 
s.t.  <  ®i(yu)  <  o,  Vi  =  i  •  •  •  |i| 
l^-(y«)  =  o,Vj  =  i---|E| 


(4) 


Eq  (4)  is  no  longer  a  Rayleigh  quotient,  solvable  by  general  eigen-decomposition,  as  compared 
to  the  traditional  unconstrained  NCUT  formulation.  As  such,  there  is  no  closed  form  solution  for  this 
problem.  In  particular,  [6]  showed  that  it  is  NP-hard,  in  general,  to  solve  for  general  eigenvectors 
under  linear  inequalities.  So,  we  propose  an  iterative  algorithm  to  find  the  global  minimum  of  this 
optimization  problem  using  Dinkelbach’s  method  for  fractional  programming  (FP)  [8,21].  To  make 
the  paper  self-contained,  we  give  a  brief  description  of  this  algorithm  next. 


4.1  Dinkeibach  Algorithm  for  Fractional  Programming 

Given  two  continuous  functions  /  :  R71  — >  R  and  g  :  — ►  R  defined  on  a  compact  set  S  C  such 

that  g  (x)  >  0  Vx  E  S ,  the  fractional  programming  problem  is  to  find  the  global  minimizer,  x*,  of 
h  (x)  =  According  to  the  parametric  approach  of  Dinkeibach  [8],  x*  minimizes  this  problem 
if  and  only  if  JT(x*,  A*)  =  min  [f  (x)  —  A *g  (x)],  where  A*  =  h  (x*).  Dinkeibach  proved  that 
jr(x,  A)  is  monotonically  decreasing  in  A.  This  equivalence  was  extended  to  prove  that  A*  can  be 
reached  iteratively.  In  fact,  he  proposed  an  algorithm  that  produces  a  sequence  of  monotonically 
decreasing  values  of  A^  =  h  (xW),  which  converges  superlinearly  to  A*. 

The  Dinkeibach  algorithm  was  extended  by  [21]  to  provide  a  general  framework  for  FP’s,  sum¬ 
marized  below  as  Algorithm  1.  A*  is  the  global  minimum  value  of  the  objective  function.  Here, 
we  emphasize  that  the  superlinear  convergence  property  is  only  regarding  the  iterations  needed  to 
achieve  A*  and  not  the  convergence  of  each  iteration.  Note  that  each  iteration  involves  solving  a 
different  optimization  problem,  which  might  be  an  NP  hard  problem  in  its  own  right!  Consequently, 
choosing  the  initial  guess  x(°)  is  not  a  trivial  task,  since  setting  it  to  an  arbitrary  value  might  lead  to  a 
sequence  of  NP  hard  problems.  However,  we  choose  x(°)  to  reduce  the  computational  complexity  of 
each  iteration  and  the  total  number  of  iterations  required. 


4.2  Applying  the  Dinkeibach  Algorithm  to  Eq  (4) 


Eq  (4)  fits  the  form  required  to  apply  the  Dinkeibach  algorithm,  where  S  =  {yu  :  y^Rl  +  q  = 
o,  $i(y«)  <  o  Vi  =  1,  •  •  •  ,  |I|,  Vj(yu)  =  o,  Vj  =  1,  •  •  •  ,  |E|},  f  =  F,g  =  G,  and  yi0)  =  x®.  The 


resulting  optimization  problem  to  be  solved  in  each  iteration  of  Dinkelbach(S ,  F,  G,  fu\0  is  a  QP 
subject  to  convex  constraints.  Eq  (5)  shows  the  problem  to  be  solved  at  iteration  i. 
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Algorithm  1:  Dinkelbach _ 

Input  :  S,  f,  g,  x(0)  <E  <S,  e 
Output:  Ne,  {X^Q,  x*,  and  A* 

1  begin 

2  Initialization:  A(0)  =  h  (x(°))  =  ;  i  =  0;  <5(i)  <—  oo;  Ne  =  1;  A*  =  A^0^ 

3  while  S(i)  >  6  do 

4  Solve  x*  =  arg  [/  (x)  —  A ^g  (x)] 

5  <5(*  +  l)  =  |/(x*)-A«p(x*)| 

6  Ne  < —  Ne  T  1,  i  < —  i  T  1 

7  If  5{i  +  1)  <  e:  A*  =  h  (x*),  break 

8  else  A(*+1)  =  h  (x*) 

9  end 

10  end 

11  end 


y u  =  arg  min 


fu  (Q  “  A^R)  yu  +  mTyu  + 


s.t.  yueS 


(5) 


The  computational  cost  of  each  iteration  is  highly  dependent  on  the  nature  of  the  matrix  = 
Q  -  AWR.  If  MW  b  0,  the  problem  becomes  a  convex  QP,  whose  global  minimum  can  be  found 
efficiently  depending  on  the  nature  of  S.  However,  if  MW  has  at  least  one  negative  eigenvalue,  then 
finding  the  global  minimum  of  Eq  (5)  becomes  NP  hard  [19].  Therefore,  it  is  essential  that  we  study 
the  existence/uniqueness  of  an  initial  guess  yi°^  that  guarantees  the  convexity  of  each  Dinkelbach 
iteration  or  equivalently  the  positive  semi-definiteness  of  each 


4.2.1  Dinkelbach  Initialization  for  Eq  (4):  A^0) 

In  what  follows,  we  will  determine  an  a  bound  on  that  ensures  MW  t  0  V  i  =  0,  •  •  •  ,  Ne, 
thus,  making  Eq  (5)  convex  for  every  iteration  (refer  to  Theorem  1).  Furthermore,  we  show  how  to 
construct  a  valid  yi°^  that  satisfies  this  a  bound. 

Theorem  1  (a  Bound  on  A<°)).  <  a,  M«  y  0  V  i  =  0,  •  •  •  ,  iVe.  Here,  a  =  min  (^"(r^~p") 

is  a  non-trivial,  upper  bound  computed  without  eigen-decomposition. 

Proof.  We  propose  to  select  A^  in  order  to  guarantee  y  0  V  i  =  0,  •  •  •  ,  Ne.  Actually,  this  is 
equivalent  to  ensuring  >:  0,  since  there  exists  a  recursive  relationship  between  MW  and  M^: 
MW  =  M(°)  +  (A(°)  -  AW)  R.  Notice  that  R  y  0  and  A<°)  >  A®  V  i  =  1,  •  •  •  ,  Ne  (refer  to 
Section  4.1).  To  do  this,  we  study  the  relationship  between  the  eigenvalues  of  M^0),  P^,  R,  and  Q, 
which  are  matrices  in  M.NuXNu.  Here,  we  use  the  eigenvalue  notation  p &  (B)  to  denote  the  kth  largest 
eigenvalue  of  matrix  B.  Using  the  results  of  [2, 24],  we  bound  the  eigenvalues  of  M^0)  as  follows: 
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PNU  (Q)  <  Pj  (M(°))  +  A WPj  (R)  <  Pi  (Q)  V  j  =  1,  •  •  •  ,  Nu 

Since  Q,  P £>  E  ,  we  can  bound  the  eigenvalues  of  Q  in  a  similar  manner  to  give:  pNu  (Q)  > 

Pnu  (P d)  +  Pnu  (Lr)-  Here,  we  use  the  fact  that  the  smallest  eigenvalue  of  a  Laplacian  matrix  is 
zero  [23]  (i.e.  p^u  (Lw)  =  0).  This  step  is  performed  to  avoid  calculating  p^u  (Q),  which  is  compu¬ 
tationally  expensive.  But,  it  also  loosens  the  bound  on  which  may  lead  to  slower  convergence. 
However,  even  with  this  step,  our  experiments  show  that  merely  1-3  Dinkelbach  iterations  are  needed 
for  convergence. 

Combining  the  above  results,  we  find  that  pj  (M^0))  >  min  ( diag  (  Pd))  -  A(°)R(j,  j).  This  sim¬ 
plification  is  possible  because  P^  is  a  diagonal  matrix  and  its  eigenvalues  are  its  diagonal  elements 
themselves.  For  M^0)  >:  0,  we  require  that  pj  (M^0))  >0  V  j  =  1,  •  •  •  ,  Nu.  This  can  be  achieved 

when:  <  a,  where  a  =  mm  (N+Pu+N-pu) 


max  (R) 


□ 


4.2.2  Dinkelbach  Initialization  for  Eq  (4):  yi°^ 

Theorem  1  proved  the  existence  of  a  non-trivial  upper  bound  for  A*®;  however,  it  did  not  show  how 
to  find  a  particular  value  of  yi°^  that  satisfies  this  bound  (if  one  exists).  Here,  we  construct  such  a 
yi°\  by  setting  and  K  to  appropriate  values.  From  A^0)  =  H  ^yi°^  <  a,  we  obtain  a  quadratic 
feasibility  problem  with  convex  constraints,  as  shown  in  Eq  (6),  with  (Q  —  aR)  ^  0. 


We  define 


(I) :  yi0)T  (Q  -  OR)  yi0)  +  mTyi0)  +  (l-a)a-c<0 

(II) :  dTyi0)  +  <7  =  0  (d  =  Rl) 


(6) 


[(HI):  ^(yi0))  <0,^-(yD  =  0, Vi  =  !,•••  ,|I|,  Vj  =  !,■••  ,  |E| 


?(°) 


01  =  -a  \(N+)2  +  (bN_)2 
02  =  (l-  a)  [iV+  (PpJ)  +  b2N_  (lTp^ 
03  =  [(N+)2 -b(N_f 
04  =  N+  (lTp+)  -  bN _  (lTp-) 

(1  —  a)  a  —  c  =  (3\K  +  /?2 

Q  =  foK  +  Pa 


dTv(°)+  8a 

Replacing  (II)  (i.e.  K  = - ^  ®  and  enforcing  that  K  >  0,  we  obtain  an  equivalent 

feasibility  problem,  where  r  (6)  =  —  /%,  as  shown  in  Eq  (7).  Now,  all  three  constraints  are 


dependent  on  the  value  of  b  alone. 


(I): 


yi0)T  (Q  -  qR)  yi0)  +  (m  -  |^d)  yi0)  <  r  (6) 

(ii):  dTy£)+/34  <  0 

(III):  <^(yl0))  <  0,  ^(yf )  =  0,  Vi  =  1  •  •  •  |I|,  Vj  =  1 


(7) 


|E| 
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By  setting  b  to  bo  =  arg  max5>0  r  (6),  we  maximize  the  upper  bound  on  (I).  In  fact,  it  can  be  shown 

(refer  to  Appendix  A)  that  this  bound  can  be  made  arbitrarily  large,  if  bo  =  n  ,  where  n  is 

chosen  in  the  following  manner. 


K  —  1  —  £ 

Hi  =  1  H-  E 

As  stated  in  Section  2.2,  b  quantifies  how  connected  each  segment  is  to  the  rest  of  the  graph.  Ideally, 
the  value  of  b  is  dependent  on  the  final  solution  y*.  However,  subject  to  the  hard  constraints  (if  any) 
and  prior  to  computing  y* ,  bo  is  a  reasonable  estimate  of  b.  In  the  absence  of  hard  constraints  (i.e. 
N+  =  N _  =  1),  bo  — ►  1.  In  the  traditional  NCUT  formulation  [23],  no  attempts  were  made  to 
approximate  b. 

Therefore,  by  setting  b  =  bo,  we  can  construct  a  feasible  solution  to  Eq  (7)  by  solving  the  convex 
QP  shown  in  Eq  (8).  Due  to  the  convexity  of  the  cost  function  and  the  constraints  in  Eq  (8),  we  are 
guaranteed  that  a  solution  exists.  This  solution  can  be  obtained  by  using  a  suitable  QP  solver  (e.g. 
methods  based  on  trust  regions,  active  sets,  etc.)  to  find  the  global  minimum.  Moreover,  the  sparsity 
of  (Q  —  aR)  should  be  exploited  to  reduce  the  overall  computational  complexity. 


if  r  <  1 
if  r  >  1 


where 


y^0)  —  arg  min 


xT  (Q  —  aH)  x  + 


--ld 


x 


ra^x+A  <0 

S.t.  \  P 3  — 

l$i(x)  <  0,^(x)  =  0 


(8) 


4.3  Proposed  DNCUT  Algorithm 

In  Algorithm  2,  we  show  the  steps  required  to  solve  a  convexly  constrained  NCUT  problem,  under 
the  DNCUT  framework.  This  algorithm  reiterates  how  b ,  K ,  and  y^  are  computed  to  ensure  that 
>z  0  in  every  iteration  of  the  Dinkelbach  algorithm.  Notice  that  the  global  minima  of  Eq  (5)  and 
Eq  (8)  are  computed  using  the  same  convex  QP  solver.  We  gain  significant  speed  up  by  initializing 
the  solution  of  each  Dinkelbach  iteration  with  the  solution  of  the  previous  one.  In  the  discretization 
step,  we  cluster  the  values  of  the  relaxed  solution  y*  to  obtain  the  binary  solution.  Any  clustering 
algorithm  can  be  used  here  (e.g.  k-means  clustering  with  k  =  2). 

Beyond  Binary  Segmentation:  DNCUT  can  be  extended  to  multi-class  ( C  >  3)  segmentation  in  a 
recursive  fashion,  similar  to  what  was  done  in  [16,23].  We  first  construct  an  over- segmentation  of  the 
image  by  clustering  the  values  of  y*  into  k^>  C  clusters.  No  particular  choice  of  clustering  algorithm 
is  required  here.  In  our  experiments,  we  used  both  k-means  clustering  and  mean  shift  clustering  [10]. 
Then,  the  k  segments  are  greedily  merged  until  only  C  segments  remain.  At  each  step,  two  clusters 
are  merged,  if  the  NCUT  cost  between  them  is  the  largest  among  all  other  pairs  of  clusters.  This 
guarantees  that  the  merged  segments  are  the  “most  similar”.  Each  resulting  segment  does  not  need 
to  be  spatially  connected  (i.e.  spatially  fragmented  labels  might  occur).  In  our  experiments,  this 
simple  merging  method  produced  meaningful  results;  however,  more  elaborate  merging  schemes  can 
be  employed. 
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Algorithm  2:  DNCUT 


Input  :  S,  Wu,  p+,  pu ,  N+,  N-,  e,  and  e 

Output:  y*  and  A* 

t  begin 

■  =  hn  =  K.  ( £t\2 


Use  e  to  compute  b  =  bo  —  n,  j 

Store  P^,  Q,  R,  m,  and  d 
Compute  a,  (3i,  fc,  fe,  and  /?4 

Solve  Eq  (8) 

_ 


yi0)  [QP  Solver] 


Compute  K  =  —  u  y^3 ~I~P4 ,  a,  c,  q,  and  =  H  (y«°^ 


8 

9  end 


Ne,  {Aw}^f0,  y* ,  A*  j  =  Dinkelbach(S,  F  (yu),  G  (yu),  yi°},  e)  [QP  Solver] 
Discretise:  y*  = 


-biT  y*T 


4.4  Special  Case:  DNCUT  under  Linear  Constraints 

In  this  section,  we  describe  how  to  efficiently  apply  the  DNCUT  framework  under  linear  in/equality 
constraints.  The  reason  we  give  this  type  of  constraint  special  attention  is  twofold.  (1)  Linear  con¬ 
straints  encode  important  first-order  relationships  between  graph  nodes,  such  as  partial  groupings  [26] 
or  area  constraints  [9].  However,  current  methods  are  restricted  to  linear  equality  constraints.  In  what 
follows,  we  show  that  DNCUT  readily  incorporates  linear  inequalities.  (2)  There  exist  efficient  itera¬ 
tive  methods  that  solve  the  underlying  convex  QP  (e.g.  interior  point  or  active  set  methods). 

As  emphasized  earlier,  we  need  a  single  QP  solver  to  find  the  global  minima  of  the  convex  QPs 
in  Eq  (8)  and  Eq  (5).  We  use  a  basic  active  set  solver  for  the  QP  in  Eq  (9). 

min  x7  Ax  +  bTx 

/  Cx  <  d  (9) 

s.t.  < 

[Ex  =  f 

where  A  G  S^u,C  €  ]RlIlxJV“,  E  G  MdEl+1)xAr“,  b  G  RN“,d  G  and  f  G  MdEl+1).  Also,  for  our 
implementation  purposes,  we  assume  that  C  and  E  are  sparse  matrices  (e.g.  case  of  partial  groupings). 
This  assumption  highly  reduces  the  complexity  of  the  active  set  method,  whose  fundamental  step 
employs  solving  a  large,  sparse  linear  system  with  conjugate  gradients.  When  |I|  =  0,  there  exists 
a  closed  form  solution,  which  requires  solving  a  pair  of  linear  systems.  Also,  if  |E|  =  1  and  no 
hard  constraints  exist,  the  problem  degenerates  to  the  unconstrained  NCUT  problem.  When  hard 
constraints  are  applied,  the  problem  becomes  equivalent  to  interactive  NCUT.  More  details  of  the 
optimization  technique  can  be  found  in  Appendix  B. 

5  Experimental  Results 

We  conducted  a  set  of  image  segmentation  experiments  to  verify  the  correctness  of  our  formulation 
and  the  proposed  DNCUT  algorithm.  Using  the  NCUT  criterion  for  image  segmentation  was  chosen 
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to  only  validate  the  DNCUT  framework  and  to  compare  it  against  other  NCUT  algorithms.  The 
segmentation  performs  only  as  well  as  the  underlying  NCUT  formulation  (i.e.  our  objective  is  only 
a  more  general  NCUT  framework  and  algorithm).  The  results  demonstrate  that  our  algorithm  can 
accept  hard  or  other  linear  constraints  and  efficiently  derive  a  solution. 

One  problem  that  any  image-as-a-graph  approach  must  face  is  the  need  for  large  memory,  to 
accommodate  the  large  graphs  created  by  the  large  number  of  pixels  in  any  reasonably  sized  image. 
In  particular,  forming  the  weight  matrix  Wu  is  the  memory  bottle-neck  for  our  DNCUT  algorithm, 
as  it  is  for  the  original  NCUT  implementation  [23].  This  problem  may  be  partly  addressed  by  a 
multiscale,  coarse-to-fine  implementation  of  DNCUT,  similar  to  [7];  however,  its  extension  to  the 
DNCUT  framework  is  left  to  future  work.  In  our  experiments,  there  was  no  need  for  down-sampling 
images,  since  the  number  of  pixels  in  each  image  did  not  exceed  the  maximum  allowable  number  of 
nodes.  The  edge  weights  between  neighboring  nodes  are  computed  using  intervening  contours  [16], 
applied  to  the  the  grayscale  version  of  the  image.  We  use  the  default  parameters  (e.g.  scale)  that 
are  available  in  the  implementation  of  [16].  These  weights  take  on  values  ranging  from  0  to  1.  Wu 
is  made  sparse  by  setting  the  weights  between  nodes  lying  farther  than  a  certain  distance  (i.e.  10 
pixels)  from  each  other  to  zero.  For  an  unknown  node  z,  p+  (z)  and  p“  (z)  are  heuristically  set  to  the 
maximum  similarity  between  node  z  and  all  nodes  of  Sa  and  Sb  respectively.  In  the  absence  of  hard 
constraints  (i.e.  N+  =  N _  =  1),  we  set  p+(z)  =  p“(z)  =  p.  In  fact,  the  DNCUT  solution  becomes 
equivalent  to  the  traditional  NCUT  solution  when  p  — >  0.  For  our  experiments,  we  set  p  =  Here, 
we  note  that  if  a  probabilistic  model  exists  for  the  graph  nodes,  p+  (z)  and  p~  (z)  can  be  set  to  the 
likelihood  that  node  z  belongs  to  A  and  B ,  respectively.  The  tolerance  values  are:  e  =  10-3  and 

e  =  ur5. 

In  what  follows,  we  give  the  complexity  of  the  DNCUT  algorithm  (Section  5.1),  compare  the 
performance  of  DNCUT  to  two  other  NCUT  algorithms  when  applied  to  the  unconstrained  binary 
NCUT  problem  (Section  5.2),  and  show  segmentation  results  of  DNCUT  when  applied  to  the  linearly 
and  nonlinearly  constrained  binary  NCUT  problem  (Section  5.4)  and  to  the  unconstrained  multi-class 
segmentation  problem  (Section  5.5). 

5.1  Complexity: 

All  our  image  segmentation  experiments  were  executed  using  MATLAB  7.6  on  a  2.4  GHz,  4GB  RAM 
PC.  With  hard  and  linear  constraints,  our  algorithm  required  1-3  Dinkelbach  iterations  to  converge. 
In  general,  the  worst  case  complexity  of  our  algorithm  is  O  where  Nu  is  the  total 

number  of  unknown  nodes  in  Q  and  O  (Nu2^  is  the  worst  case  complexity  of  solving  a  sparse  linear 
system  with  Nu  variables  using  conjugate  gradients,  po  is  the  maximum  number  of  Dinkelbach 
iterations  required  for  convergence  and  pa  is  the  maximum  number  of  active  set  iterations  needed  for 
one  Dinkelbach  iteration  to  converge.  When  the  constraints  are  linear  equalities,  pa  =  1. 

5.2  Validation: 

Here,  we  demonstrate  the  correctness  of  DNCUT  by  comparing  it  to  two  previous  NCUT  implemen¬ 
tations  [7, 23]  applied  to  the  same  image.  This  is  done  quantitatively  for  the  case  of  unconstrained 
NCUT  (i.e.  |E|  =  |I|  =  0),  where  the  global  solution  to  the  relaxed  NCUT  problem  is  known  to  be 
the  second  smallest  generalized  eigenvector  (i.e.  eig2(D  —  W,D)).  In  this  case,  N+  =  N-  =  1 
(i.e.  “source”  and  “sink”  are  only  included).  For  all  three  algorithms,  we  explicitly  apply  the  hard 


H(0) 
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constraints  by  setting  y+  =  1  and  y_  =  —  V 

The  NCUT  methods  we  compare  against  are  the  original  algorithm  described  in  [23]  and  its  multi¬ 
scale  version  described  in  [7].  These  algorithms  use  implicitly  restarted  Arnoldi/Lanczos  methods  for 
sparse  matrices  to  perform  eigen-decomposition  in  an  iterative  fashion  [13].  The  run-time  of  this  it¬ 
erative  decomposition  is  determined  by  the  relative  tolerance  eeig.  The  iterations  are  terminated  when 
the  relative  change  in  eigen  solutions  at  the  current  and  previous  iterations  is  less  than  eeig.  On  the 
other  hand,  DNCUT  solves  the  unconstrained  problem  by  solving  a  pair  of  sparse  linear  systems.  In 
fact,  this  pair  of  systems  can  be  solved  simultaneously,  using  the  iterative  conjugate  gradient  method. 
This  is  true  since  the  conjugate  directions  for  one  of  the  systems  can  be  used  for  the  other.  The  run¬ 
time  of  DNCUT  is  determined  by  the  relative  tolerance  ccg?  which  defines  the  stopping  criterion  for 
the  conjugate  gradient  method.  When  the  relative  change  in  the  solution  to  the  linear  system  is  less 
than  6cg>  the  algorithm  terminates.  For  all  three  algorithms,  the  run-times  do  not  include  the  time 
needed  to  construct  the  weight  matrix  Wu. 

In  Figure  5,  we  show  comparative  results  for  the  three  aforementioned  algorithms  when  applied 
to  images  in  the  Berkeley  segmentation  dataset  [17].  Each  image  in  this  dataset  has  ~155,  000  pix¬ 
els/nodes.  We  aim  to  show  how  the  NCUT  cost  (as  defined  in  Eq  (4))  varies  for  each  algorithm  and  for 
different  run-times  (i.e.  for  different  stopping  criteria).  On  each  image  in  the  dataset,  we  ran  the  three 
algorithms  with  varying  stopping  criteria  and  registered  their  corresponding  NCUT  costs.  H  (O), 
H  (M),  and  H  (D)  denote  the  costs  of  the  original  NCUT  [23],  multiscale  NCUT  [7],  and  DNCUT 
algorithms  respectively.  A  total  of  20  stopping  criteria  were  used  for  each  algorithm.  Figure  5  plots 
the  three  NCUT  costs  at  each  run-time  t.  This  cost  is  averaged  over  all  images  in  the  dataset.  We  also 
show  the  standard  deviations  of  these  measurements.  Here,  we  used  linear  interpolation  to  complete 
the  plots.  It  is  obvious  that  as  the  run-time  of  an  algorithm  increases,  its  NCUT  cost  decreases  till 
it  reaches  a  stable  value.  All  three  algorithms  exhibit  this  variation.  From  the  plots,  we  conclude 
that  the  original  algorithm  outperforms  the  multiscale  one  by  10.9  dB ,  while  the  DNCUT  algorithm 
outperforms  the  original  one  by  3.7  dB.  This  points  to  the  obvious  fact  that  solving  Eq  (4)  directly 
(i.e.  NCUT  on  the  augmented  graph)  is  not  equivalent  to  solving  the  NCUT  problem  on  the  pixel 
nodes  alone.  Moreover,  as  p  — »  0,  the  original  NCUT  solution  will  approximate  the  DNCUT  one 
and  H  (D)  — ►  H  ( O ).  Ideally,  the  comparison  should  not  be  a  relative  comparison  between  the  three 
methods,  but  instead  a  comparison  of  each  method  with  the  global  minimum  of  the  original  NCUT 
problem  in  Eq  (2).  Since  this  problem  is  NP  hard,  obtaining  its  global  minimum  for  non-trivially 
sized  problems  is  infeasible. 


time  (seconds)  lime  (seconds)  time  (seconds) 


Figure  5:  NCUT  costs  for  the  three  NCUT  algorithms  are  plotted  versus  run-time.  The  unconstrained 
NCUT  problem  is  addressed  here.  The  red  values  are  averaged  over  all  the  images  in  the  Berkeley 
segmentation  dataset.  The  standard  deviations  are  also  plotted  as  blue  bars. 
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We  also  provide  four  qualitative  examples  in  Figure  6.  Columns  (b)-(d)  show  the  relaxed  solutions 
(i.e.  prior  to  any  discretization)  produced  by  the  original,  the  multiscale,  and  the  DNCUT  algorithms 
on  sample  images  in  column  (a),  respectively.  These  solutions  were  obtained  after  the  stable  NCUT 
value  for  each  algorithm  was  reached.  Comparing  the  segmentation  results,  we  see  that  the  DNCUT 
solution  is  more  detailed,  which  facilitates  segmentation.  If  we  consider  the  bird  image  in  the  second 
row,  we  see  that  the  DNCUT  solution  plays  the  role  of  a  soft  labeling,  where  pixels  with  similar 
values  are  grouped  together.  So,  the  sky  pixels  are  much  darker  than  the  foreground.  This  is  not  the 
case  for  the  other  two  algorithms,  as  they  do  not  utilize  the  augmented  graph. 


(a):  original  image 


(b):  original  NCUT  (c):  multiscale  NCUT 


(d):  DNCUT 


Figure  6:  columns  (b)-(d)  show  the  stable  NCUT  solutions  yielded  by  the  three  methods  respectively, 
when  applied  to  the  images  in  column  (a).  The  unconstrained  NCUT  problem  is  addressed  here. 


5.3  Computational  Analysis: 

Here,  we  focus  on  a  computational  analysis  of  the  three  NCUT  methods,  when  applied  to  the  prob¬ 
lem  of  unconstrained  NCUT  for  the  images  in  the  Berkeley  dataset.  We  study  the  relative  change  of 
the  NCUT  solution  with  run-time.  At  every  run-time  t,  we  calculate  the  relative  solution  change  as: 
Ae  =  .  This  change  is  averaged  over  all  the  images  in  the  dataset.  In  Figure  7,  we  plot 

Ae,  as  a  percentage,  for  each  algorithm.  All  three  algorithms  show  the  same  type  of  variation.  The 
initial  monotonic  increase  in  Ae  is  followed  by  a  monotonic  decrease  till  a  stable  solution  is  achieved. 
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Original  NCUT  Multiscale  NCIIT  DNCUT 


Figure  7:  Relative  change  of  the  NCUT  solution  with  run-time 


Moreover,  DNCUT  converges  to  a  stable  solution  much  quicker  than  the  other  two  algorithms.  For 
a  5%  change,  the  original  NCUT  algorithm  requires  65  seconds  on  average  to  stabilize,  while  the 
multiscale  algorithm  requires  about  35  seconds.  On  the  other  hand,  DNCUT  requires  about  11  sec¬ 
onds  to  reach  a  stable  solution.  The  disparity  between  the  three  algorithms  is  due  to  the  inherent 
computational  nature  of  each  algorithm.  From  these  empirical  results,  we  conclude  that  DNCUT  (or 
equivalently  the  conjugate  gradient  algorithm)  has  significantly  better  convergence/stability  proper¬ 
ties  than  the  original  or  multiscale  NCUT  algorithms  (or  equivalently  the  Arnoldi/Lanczos  method 
for  eigen-decomposition). 


5.4  Interactive,  Linearly  Constrained,  and  Nonlinearly  Constrained  NCUT: 


Here,  we  conducted  three  experiments,  whereby  a  different  type  of  constraint  on  the  nodes  of  the  aug¬ 
mented  graph  is  applied  in  each  experiment.  The  first  experiment  addresses  the  problem  of  interactive 
NCUT.  In  Figure  8,  we  consider  the  case  of  interactive  NCUT,  where  a  user  marks  the  hard  constraints 
on  the  displayed  image  in  column  (b)  with  green  strokes  defining  Sa  and  red  strokes  defining  Sb- 
Columns  (c)  and  (d)  show  the  segmentation  results  produced  by  DNCUT,  with  and  without  the  hard 
constraints  respectively.  Column  (e)  shows  the  corresponding  segmentations  produced  by  the  original 
NCUT  algorithm  without  the  hard  constraints  (i.e.  unconstrained  binary  segmentation).  It  is  evident 
that  with  additional  (user-defined)  constraints,  the  resulting  binary  (i.e  foreground  vs.  background) 
segmentations  become  more  perceptually  relevant. 

The  second  experiment  addresses  the  problem  of  linearly  constrained  NCUT  and  shows  some 
results  in  Figure  9.  In  this  experiment,  we  apply  two  types  of  linear  constraints:  (1)  partial  pixel 
groupings  and  (2)  —  6q1  <  y u  <  1  to  produce  binary  segmentations  shown  in  column  (c).  The  box 
constraint  (2)  is  a  linear  relaxation  of  the  original  discrete  constraint  in  Eq  (2).  Such  linear  inequalities 
cannot  be  handled  by  other  NCUT  algorithms.  Unconstrained  binary  segmentations  produced  by 
DNCUT  and  the  original  NCUT  algorithm  are  presented  in  columns  (d)  and  (e)  respectively.  The 
partial  groupings  are  determined  by  user  defined  strokes  in  column  (b).  The  same  colored  pixels 
define  nodes  of  the  graph  that  should  belong  to  the  same  segment. 

As  in  the  second  experiment,  the  third  one  also  addresses  the  problem  of  partial  groupings  on 
the  nodes  of  the  augmented  graph.  However,  in  this  version  of  the  partial  grouping  problem,  the 
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box  constraint  on  yu  is  replaced  by  a  ball  constraint  (i.e. 

This  ball  constraint  is  a  relaxation  on  the  box  constraint,  since  it  does  not  imp 
yu(z)  E  [—bo,  l]Vz  =  Nu.  This  is  an  example  of  how  a  nonlinear  (quadratic)  convex  con- 
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straint  can  be  applied  to  the  DNCUT  framework.  We  use  an  interior  point  (barrier)  method  to  solve 
the  DNCUT  problem  under  these  two  constraints.  Since  there  is  a  single  nonlinear  constraint,  solving 
this  version  of  the  partial  grouping  problem  is  more  efficient  than  the  version  in  the  previous  experi¬ 
ment.  Figure  10  shows  some  segmentation  results  on  real  images.  The  constrained  DNCUT  solutions 
(after  discretization)  are  shown  in  column  (c),  while  the  unconstrained  binary  segmentations  produced 
by  DNCUT  and  the  original  NCUT  algorithm  are  presented  in  columns  (d)  and  (e)  for  comparison. 
The  partial  groupings  are  determined  by  user  defined  strokes  in  column  (b).  The  same  colored  pixels 
define  nodes  of  the  graph  that  should  belong  to  the  same  segment. 


(a):  original  image 


(b):  hard  constraints 


(c):  interactive  NCUT  (d):  uncon.  DNCUT 


(e):  uncon.  NCUT 


Figure  8:  Shows  examples  of  interactive  NCUT,  as  compared  to  unconstrained  NCUT  and  unconstrained  DNCUT. 
The  original  images  are  shown  in  column  (a).  In  Column  (b),  we  show  the  hard  constraints  of  Sa  (green)  and  Sb 
(red),  as  marked  by  a  user.  Using  these  hard  constraints,  column  (c)  shows  the  interactive  segmentations  produced 
by  DNCUT.  Columns  (d)  and  (e)  display  the  binary,  unconstrained  NCUT  segmentations,  produced  by  DNCUT  and 
the  original  NCUT  algorithm  respectively.  These  last  two  columns  are  shown  for  comparison  with  the  interactive 
DNCUT  results. 


5.5  Multi-Class  Segmentation: 

Figure  11  shows  examples  of  unconstrained  multi-class  segmentation  where  C  =  2, 3, 4.  (b)  shows 
the  DNCUT  solution  to  the  unconstrained  problem.  The  clustering  and  merging  algorithm  in  Section 
4.3  is  used  to  produce  the  segmentations  in  (d),  (f),  and  (h),  where  the  corresponding  boundaries  are 
drawn  in  (c),  (e)  and  (g).  Even  though  this  multi-class  segmentation  algorithm  is  suboptimal,  it  results 
in  reasonable  segmentations.  However,  its  performance  is  correlated  with  the  quality  /detail  of  pixel 
groupings  in  the  DNCUT  solution.  For  example,  in  the  FLOWERS  case,  the  flowers  are  significantly 
delineated  in  the  DNCUT  solution;  however,  the  lady  bugs  on  the  flower  petals  and  the  blades  of 
grass  in  the  out-of-focus  background  are  not.  This  is  primarily  due  to  the  nature  of  the  weight  matrix 
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(a):  original  image  (b):  hard  constraints  (c):  con.  DNCUT  (d):  uncon.  DNCUT  (e):  uncon.  NCUT 


Figure  9:  Shows  examples  of  linearly  constrained  NCUT,  as  compared  to  unconstrained  NCUT  and  unconstrained 
DNCUT.  Two  types  of  linear  constraints  are  applied:  partial  pixel  groupings  and  a  box  constraint  on  the  values  of  the 
DNCUT  solution  (i.e.  —  b^l  <  yn  <  1).  Partial  groupings  are  marked  by  users  in  red  and  green  strokes,  as  shown  in 
column  (b).  Column  (c)  shows  the  binary  DNCUT  solutions  to  the  linearly  constrained  NCUT  problem.  For  visual 
comparison,  columns  (d)  and  (e)  show  the  binary,  unconstrained  NCUT  segmentations,  produced  by  DNCUT  and 
the  original  NCUT  algorithm  respectively. 


used.  In  fact,  our  experiments  used  intervening  contours  at  a  single  scale  and  no  color  information 
was  exploited.  Incorporating  more  visual  cues  (e.g.  color)  into  the  edge  weights  and  finding  efficient 
ways  to  combine  results  at  multiple  scales  are  left  for  future  work. 

Next,  we  show  how  the  low-level,  multi-class  segmentation  results  for  the  three  algorithms  relate 
to  human  segmentation  results.  Here,  we  use  the  Berkeley  segmentation  dataset,  which  contains 
multiple  benchmark  (human)  segmentations  for  each  image  in  the  dataset.  To  each  image,  we  applied 
the  three  multi-class  algorithms,  where  the  number  of  classes  was  set  to  the  number  of  segments 
labeled  in  the  benchmark  segmentations  corresponding  to  this  image.  Then,  the  resulting  (binary) 
segment  boundaries  are  averaged  over  all  the  benchmark  segmentations.  We  used  the  evaluation 
kit  of  [17]  to  plot  the  precision-recall  curves  shown  in  Figure  12.  The  three  algorithms  yield  very 
similar  F-scores,  which  are  higher  than  random  performance  (0.41)  and  significantly  less  than  human 
performance  (0.79).  They  rank  low  on  the  list  of  state-of-the-art  segmentation  algorithms.  More 
importantly,  comparing  these  NCUT  algorithms  together,  we  see  that  the  multiscale  one  has  the  worst 
F-score  (0.46),  while  the  original  one  has  the  best  F-score  (0.54).  DNCUT  registers  an  F-score  (0.52) 
that  is  very  similar  to  the  original  algorithm.  This  is  the  case,  even  though  DNCUT  was  shown 
(in  Section  5.2)  to  yield  a  smaller  NCUT  cost  than  the  other  algorithms.  This  discrepancy  points 
to  the  fact  that  the  NCUT  formulation  for  image  segmentation  does  not  correlate  well  with  human 
segmentation. 
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(a):  original  image 


(b):  constraints 


(c):  con.  DNCUT 


(d):  uncon.  DNCUT 


(e):  uncon.  NCUT 


Figure  10:  Shows  examples  of  nonlinearly  constrained  NCUT,  as  compared  to  unconstrained  NCUT  and  uncon¬ 
strained  DNCUT.  Two  types  of  constraints  are  applied:  partial  pixel  groupings  and  a  ball  constraint  on  the  values 


of  the  DNCUT  solution  (i.e. 


Yu 


M 


<  nh 


1+bp 

2 


).  Partial  groupings  are  marked  by  users  in  red  and 


green  strokes,  as  shown  in  column  (b).  Column  (c)  shows  the  binary  DNCUT  solutions  to  the  nonlinearly  constrained 
NCUT  problem.  For  visual  comparison,  columns  (d)  and  (e)  show  the  binary,  unconstrained  NCUT  segmentations, 
produced  by  DNCUT  and  the  original  NCUT  algorithm  respectively. 


6  Conclusions  and  Future  Work 

We  have  presented  a  unifying  DNCUT  framework  for  solving  convexly  constrained  NCUT  problems 
with  data  priors  on  the  augmented  graph  nodes.  We  avoid  using  traditional  eigen-decomposition,  due 
to  its  restrictions  on  the  types  of  constraints  it  can  handle  and  its  computational  complexity.  In  this 
framework,  any  convexly  constrained  NCUT  problem  can  be  converted  into  a  sequence  of  convex 
QP’s.  In  the  case  of  linear  constraints,  we  propose  an  algorithm  to  efficiently  find  the  global  so¬ 
lution  of  each  QP.  To  validate  the  correctness  of  DNCUT,  we  compare  it  to  state-of-the-art  NCUT 
algorithms.  As  compared  to  these  algorithms,  DNCUT  provides  a  better  and  more  computationally 
efficient  solution.  We  also  show  results  of  binary  segmentation  under  hard  and  linear  constraints,  in 
addition  to  multi-class  segmentation  results.  In  the  future,  we  plan  to  develop  a  multiscale  version 
of  this  framework  to  handle  larger  graphs  and  to  incorporate  grouping  information  from  weight  ma¬ 
trices  computed  at  different  scales.  Furthermore,  we  plan  to  improve  the  multi-class  segmentation 
extension. 
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TIGER 


(a):  image  (b):  DNCUT 


(c):  C  =  2  (d):  C  =  2 


(e):  C  =  3  (f):  C  =  3 


(g):  C  =  4  (h):  C  =  4 


OWLS 


(a):  image  (b):  DNCUT 


(c):  C  =  2  (d):  C  =  2 


(e):  C  =  3  (f):  C  =  3 


(g):  C  =  4  (h):  C  =  4 


FLOWERS 


(a):  image  (b):  DNCUT 


(c):  C  =  2  (d):  C  =  2 


(e):  C  =  3  (f):  C  =  3 


(g):  C  =  4  (h):  C  =  4 


Figure  11:  Examples  of  unconstrained  multi-class  segmentation  with  C  —  2, 3, 4,  when  applied  to  three  images: 
TIGER ,  OWLS ,  and  FLOWERS,  (b)  shows  the  DNCUT  solution  to  the  unconstrained  problem,  (c),  (e),  and  (g) 
show  the  boundaries  of  the  labeled  segments  in  (d),  (f),  and  (h)  respectively.  We  refer  the  readers  to  [1]  for  all  the 
segmentation  results. 


Figure  12:  Precision-Recall  curves  for  the  three  NCUT  algorithms,  when  applied  to  the  Berkeley 
human  segmentation  dataset. 
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APPENDIX  A 

Dinkelbach  Initialization  for  Eq  (4)  [Section  (4.2)] 

In  this  part,  we  will  give  a  more  detailed  description  of  bo  =  arg  maxfe>0  r  ( b ).  We  expand  r  ( b )  as  in 

N-  ( 

Eq  (10),  where  r  = 


N+  iTvl 


r  (b)  = 


—a 

(iV+)2  +  (WV_)2 

N+  | 

1  -  bN_  | 

f 

(iV+)2  —  b  (N-)2 

-(!-«)  N+  ( 1TPJ  )  +  bzN-  ( 1 


rr 


ab  [b  +  1] 


Cp  ujN*. 

nz 


[l-r] 


>  -b 


N, 


(rT^+)+62iv_(iT^)' 


(10) 


Since  a,  b ,  lTpu ,  N+,  N- ,  and  r  are  non-negative  numbers,  we  deduce  that  r  ( b )  can  be  arbitrarily 
large  by  setting  6  =  =  ,  where  k  is  chosen  such  that, 


k,  =  1  —  e  if  r  <  1 

where 

K  =  l  +  £  if  T  >  1  |  T  — 


£  >  0,  £  ^  0 

N-  ( lTp + 


N+  lTp~ 


APPENDIX  B 

Solving  a  Convex  QP  with  Linear  Constraints  [Section  (4.4)] 

In  the  next  three  parts,  we  will  describe  how  to  solve  a  convex  quadratic  programming  problem  with 
general  linear  constraints,  of  the  following  form: 


min  x7  Ax  +  b7  x 

fCx  <  d 
s.t.  <  ~  _ 

I  Ex  =  f 


(11) 


Next,  we  consider  the  case  when  only  a  single  equality  constraint  exists  (e.g.  in  the  case  of  uncon¬ 
strained  or  interactive  NCUT).  Then,  we  consider  the  case  of  multiple  equality  constraints  followed 
by  general  linear  constraints  (i.e.  in/equalities). 
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Solving  a  Convex  QP  with  a  Single  Linear  Equality  Constraint 

Here,  we  will  provide  the  closed  form  solution  (x*)  to  the  following  convex  QP  problem. 


min 


xT  Ax  +  bTx 


s.  t.  eTx  +  /  =  0 


(12) 


We  derive  the  Lagrangian  dual  of  Eq  (12),  £  (x,  u),  and  determine  the  primal-dual  solution,  x*  and 
z/\  in  closed  form.  Equations  (14,15)  evaluate  these  optimal  solutions.  Note,  that  the  closed  form 
solution  of  the  primal  problem  is  valid,  since  C  (x,  v)  is  convex  in  x.  Also,  the  form  of  i/*  directly 
follows  from  the  concavity  of  £(x*,z/).  The  optimal  primal  and  dual  solutions  can  be  obtained 
efficiently  by  solving  a  pair  of  linear  systems:  Ax  =  b  and  Ax  =  e,  using  preconditioned  conjugate 
gradients. 


£  (x,  v) 

Primal  Solution:  x* 

Dual  Solution:  v* 


x  +  uf 


arg  min  £  (x,  v) 

X 


b  +  u*e 


arg  max  £  (x* ,  v) 

V 


2/  -  erA  1b 
eTA_1e 


Solving  a  Convex  QP  with  Multiple  Linear  Equality  Constraints 

Here,  we  will  provide  the  closed  form  solution  (x*)  to  the  following  convex  QP  problem. 


(13) 

(14) 

(15) 


xT  Ax  +  bTx 


mm 

s.  t.  Ex  +  f  =  0 


Similar  to  the  previous  part,  we  solve  the  primal  dual  problems  by  solving  the  linear  system  in  Eq 
(16).  However,  in  this  case,  the  primal  and  dual  variables  are  coupled  in  the  same  linear  system. 


'  2A 

Et  ' 

'  x*  ' 

b 

E 

0 

V* 

f 

(16) 


Solving  a  Convex  QP  with  General  Linear  Constraints 

Here,  we  will  discuss  an  active  set  based  method  that  iteratively  solves  the  following  convex  quadratic 
programming  problem.  Let  /  be  the  set  of  inequality  constraints. 


min  x7  Ax  +  b7  x 

fCx  <  d 
s.t.  <  -  _ 

|  Ex  =  f 


We  present  the  main  steps  involved  in  applying  the  active  set  method. 
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1.  Initialize  a  counter  k  —  0.  Choose  a  feasible  initial  guess  x®,  which  satisfies  a  certain  set 

(k) 

SE  j  of  the  constraints  with  equality,  including  all  the  equality  constraints  and  a  subset  of  the 
inequality  constraints.  These  constraints  are  referred  to  as  active  constraints.  They  correspond 
to  the  following  set  of  linear  equalitiesiG^x  =  h^k\ 

2.  Compute  a  step  direction  5^k\  by  solving  a  convex  QP  with  linear  equality  constraints  corre¬ 
sponding  to  those  of  SE^ .  This  is  done  using  the  method  outlined  in  the  previous  part. 

5^  =  arg  min  zTAz  +  +  2Ax^^  z 

s.t.  G^z  =  0 

3.  If  5^  =  0,  then  check  the  Lagrangian  multipliers  corresponding  to  the  active  inequality  con¬ 

straints.  If  they  are  all  positive,  then  the  final  solution  has  been  obtained;  otherwise,  remove  the 
constraints  corresponding  to  the  negative  multipliers  from  S .  On  the  other  hand,  if  5^  /  0, 
take  a  step  from  the  previous  solution  along  5^  (i.e.  —  x(fc+1)  +  The  step  size 

is  computed  as  follows,  where  is  the  zth  row  of  G®  and  h(i)W  is  the  ith  element  of  h^k\ 


a(k)  _  mjn 


gf)Tx(fc) 


ik)T5(k) 


where  p  —  arg  min 

i  e  S(k)  n  I 

gf  )t5^  >  0 


—  gf)Tx(fc) 
g pr5{k) 


This  step  value  guarantees  that  the  current  solution  lies  on  the  set  of  active  constraints.  If 
<  1  ,  then  add  the  pth  constraint  to  SE\  Increment  k. 

4.  Iterate  through  steps  2-3  until  convergence. 
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